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We investigate the exact solution of BCS pairing model using direct diagonalization of Fock 
space. By the data analysis and numerical calculation, we verify the symmetry between energy 
spectrum of Fock subspaces, obtain the common structure features of energy gaps and energy bands 
in Hamiltonian spectrum of reduced model, propose the formula to estimate the lowest energy levels 
in all of the subspaces of reduced model, and suggest a scheme to estimate the respective energy 
spectrum which can reveal the structure of energy spectrum of pairing model. 

PACS numbers: 74.20.Fg, 03.65.-w, 03.67.Lx 

Recently, the Bardeen, Cooper and Schrieffer (BCS) model for superconductivity pj has been connected with the 
problems in different areas of physics such as superconductivity, nuclear physics, physics of ultrasmall metallic grains 
and color superconductivity in quantum chromodynamics 0-L| Since Ricardson's works from 60's to now , the 
exact solution of the reduced (constant) BCS pairing model has been well known. However, for the more general cases 
in which the coupling coefficients among different pairs are different, the exact solution of BCS pairing model is on 
studying I* 1 fact, the exactly solvable models have been proven to be very useful tools to understanding 

the physics of strongly correlated many-body quantum systems. There are a serial of the important and interesting 
works about the pairing model being published that have not cited here (see the references in 0,0). 
Let us consider the BCS pairing Hamiltonian Q: 
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where n± m = c* ±m c± m are the electron number operators, c}„(c m ) is the fermionic creation (annihilation) operator, 
and the coupling coefficients are real and symmetric, that is V m i = Vi m . Based on the Refs. 0, 01 111- one is able to 
study equivalently its spin-analogy form as the following 
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where e m = [e m — Bp) — V mm , and the constant term X^m 6 ™/^ has been ignored, which vanishes anyway since we 
cut off symmetrically above and below ef- Note that L is the number of the pairs or qubits here. 

Quantum theory tell us, it is natural to exactly solve the pairing model by the direct diagonalization in Fock 
space. Several works have been published in this aspect, for example 0, |nj. However, its feasibility still has to 
be further considered under no any approximation. Clear calculations and interesting analysis of the results need 
to further study. Obviously, a difficulty of exact solution of spin-analogy of pairing model is, in the large-L limit, 
the obsession of exponentially complicated problem to the direct diagonalization. Fortunately, this difficulty can be 
partially avoided by making use of decomposition of Fock space. At least, for some given subspaces of Fock space, 
the direct diagonalization becomes a polynomial problem. On the other hand, the direct diagonalization needs that 
. , one knows the general forms of the Hamiltonian in the given subspaces. This implies that a key problem is to seek for 
their explicit forms. As soon as they are found and expressed, it is indeed feasible to exactly solve the pairing model 
by the direct diagonalization of Fock space. 

Actually, our motivation is arose originally by our study on quantum simulation of pairing model on a quantum 
computer (HHEI. 

Because we did not know how to implement the relevant approximations in the simulating 
solution of pairing model in a quantum computer, and we realize that the direct diagonalization will be feasible and 
direct at least on a quantum computer, and so we would like to develop it for the quantum simulation of pairing 
model in the near future. 
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Now we first recall the technology of the decomposition of Fock space and give out a strict and simple proof about 
it. In fact, the spin space S^ n in a L-pair (qubit) system can be divided into the different subspaces which correspond 



to the different numbers of spin-up states, that is 5, 
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0. Therefore, we have proved that the pairing model Hamiltonian is able to be decomposed 



into the direct sum of submatrices in Fock subspaces, that is 
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It is clear that two the simplest eigenvectors of Hp L ^ are in the one dimension Sq^' and Sj^' respectively with the 
eigenvalues =F X)m=i e m/2. Thus, the concerned subspaces or submatrices of Hamiltonian in our method only includes 
those from 1 to L — 1. 

Furthermore, we can derive out 
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where -FfJ~, 2 ...j -j - r is such a permutation that i c is dropped, r is added, and {ii,i2, ••• ,* c -i,*c+i, ■•• is 
rearranged according to their values from small to large. In terms of the orthogonality of bases, we obtain immediately 
the general and explicit forms of Hamiltonian submatrices of pairing model in the Fock space. Consequently, we can 
carry out the direct diagonalization of them within the limitation of the computer's power. 

It must point out that in the half filling case i.e. n — [L/2], the dimension of its submatrix H^ hi is the highest 
and arrives at L\/([L/2]\[(L + l)/2]!), where U [M]" means taking the integer part of M. Note the fact that H^ hi is 
a sparse matrix, for example, the dimension of is 12870, but the number of its nonzero elements is only 65 x 65, 

one is able to diagonalize it to some L. Nevertheless, with L increasing, the dimension of half filling submatrix 
will exceeds quickly the limits of large scale diagonalization in a classical computer. Perhaps, quantum simulation 
can solve this difficulty 0, E3, 

Here, we focus on the problem how obtain the more knowledge of Hamiltonian 
spectrum for the moderate L, which is helpful for understanding many body quantum theory, as well as for providing 
a precision comparison with the possible result of quantum simulating in the near future. Actually, we note that 
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Obviously, it is easy to diagonalize the ff^y for L = 100 
strictly obtain the general and explicit form of H^t,i an d the arbitrary subspace 



180 in an ordinary PC 0|. By making use of eq.Q we 

In this letter, we consider 

the typical three kinds of models, that is, the reduced (constant) model in which Vij = V (constant) for any i and 
j, the nearest neighbor model in which Vij = Su^ xV (only there is the coupling in the nearest neighbor two pairs) 
and the third model which takes Vij — (3V/\j — i\ (the coupling coefficient decreases pro rate with the "distance"). 



They are denoted respectively by the subscripts RE, NB and RA. Moreover, we take V = 2 x 10 



iV/X, 



(3 = IQ- 1 and A is a running parameter which can be taken as 10, 20 and etc in the following numerical calculation. 
Figffl shows the energy spectrum of subspace 1 of H p which is obtained by the direct diagonalization. 
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FIG. 1: The energy levels in the three kinds of models RE, NB and RA (L = 170, A = 20) 



Fig.l and the other figures in this letter are all so-called energy level figure, in which, every line segment represents 
an energy level (or several energy levels whose number is equal to its degeneration degree), that is, an eigenvalue of 
Hamiltonian, and y-axis indicates the values of energy levels. In fact, for saving space, we always put a set of energy 
level figures together, where every sub-figure occupies one column and on its top a name of the considered models 
(Fig.l) or an ordinal (not including the trivial ones of Sq and S^) of subspace (the others) are denoted. It must 
be emphasized that in our energy level figures, there are some energy levels that have not been distinguished clearly 
since the limitation of scale of y-axis and resolution from printer or monitor. Actually, it is just the effect we want to 
display, because those dense energy levels (large number of energy levels within the unit scale of energy) form so-called 
"energy bands" (the line segment with some width), two near but obviously separate energy levels (the difference of 
two near energy levels has the larger value) forms so-called "energy gaps" (two parallel line segment with an obviously 
larger interval). For example, it is easy to see that in the subspace 1 of reduced model has an energy gap standing at 
between the energy level 1 and energy level 2, and an energy band including all energy levels except for the energy 
level 1. 

By the numerical calculating and fitting, analytical continuation and then theoretical deduce, we can obtain several 
important and interesting conclusions. 

(1) Symmetry in the energy spectrum: Denoting i%E nb ra(*j n ) as the *"^ n energy level of Hp L ^ in the subspace n 
respectively for the models RE, NB and RA and arranging as E^ L \i, n) < E^ L \j, n) for i < j, our data analysis and 
numerical fitting indicate that there is the symmetry between E^ L '(i, n) and E^ L \i, L — n). That is 



(i, n) = (i, L - n) - (L + 1) - n) 
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Moreover, this symmetry is independent of the considered models here and so the subscript RE or NB or RA are 
omitted. We think that it can be an exact theoretical formula since its precision < 10 36 (16j| ■ 
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The Reduced Model (L=10} 



FIG. 2: The symmetry of energy spectrum, and the structures of energy gaps and energy bands in the reduced model, where 
A = 20. 



In fact, a p hysical symmetry should be a strict result. Newly, we have finished a strict analytical proof about this 
property (l5j|. 

(2) Structure of energy gaps and bands: From the Fig[3 it is easy to see the obvious structure of energy gaps and 
bands. Actually, we find that it is a common feature of energy spectrum in the reduced model. Here, the energy gap(s) 
means such a (some) difference(s) between two nearest neighbor energy levels that it is much larger than the other 
ones, and the energy band means such a set of energy levels with very small differences among the nearest neighbor 
energy levels, even appearing approximate degeneration. From the data of energy levels we can conclude that the 
number of energy gaps is equal to n (if n < [L/2]) or L — n (if n > [L/2]) for a given subspace n of Hp , and so is 
the number of the energy bands since it does not contain a lowest energy level in this given subspace. The positions 
of energy gaps E^(k, n) appear between E^.(L\/{{k - !)!(£ - fc + 1)!) + 1, n) and E$ {L\/((k -l)\(L-k+ l)!),n), 
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where k = 1, 2, ■ • • , n (if n < [L/2]) or k = 1,2, • •• , L — n (if n > [L/2])). We have verified, in numerical, the following 
variation rules for N < 10: the energy gaps enlarge with A or n(< [L/2]) or L increasing, and lessen with k increasing; 
the energy bands widen with k or L increasing, and narrow with A or n(< [L/2]) increasing. When n > [L/2], the 
variation of energy gaps and bands can be known by the symmetry in energy spectrum. However, FigOl shows that 
the nearest neighbor model destroys the structure of energy gaps and energy bands, and so almost does the third 
model except for the first energy gap. The destruction strength of the third model for the first energy gap depends 
on variety strength of the coupling coefficients with the distance "|j — the more rapidly the coupling coefficients 
decrease with the distance, the more largely the first energy gap is destroyed. Fig[3| gives out an example displaying 
these features. 
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FIG. 3: The symmetry of energy levels, and the destruction of the energy gaps and bands in NB and RA models, where A = 20. 

(3) Lowest energy levels: In the reduced model, we find that the lowest energy level (LEL) of Hp is at the half 
filling subspace [L /2] . In terms of the theoretical analysis on Hamiltonian and numerical fitting for the energy levels, 

we conjecture that the LEL of Hp L ^ can be estimated by the following formula with a large leading term 
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From the fact that (1, n) » (1, n - 1) + (E& ( 1 , [L/2] ) - (1, 1) )/( [L/2] - 1) - ( [L/2] - 2(n - 1) ) V obtained 



by the data analysis and numerical fitting, it follows that 
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It implies that 
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which can be used to estimate the LEL E& (1, n) in a given subspace n < [L/2]. For the other subspaces (n > [L/2]), 
we can estimate their LEL by eq.©- When L = 10, A = 20, the absolute errors of the estimated values for all of 
subspaces from 1 to [L/2] are < 10 -9 and their relative errors are < 10~ 5 . 

(4) The representative energy levels and energy spectrum: Since the typical structure of energy gaps and energy 
bands as well as the large density of energy levels for a given band in the reduced model, we can introduce a so-called 
representative energy level (REL) which is defined by an average of all of energy levels in a given energy band, that is 
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where k = 1, 2, • • • , n (if n < [L/2]) or k = 1, 2, • • • , L — n (if n > [L/2])). This REL can be used to indicate the 
corresponding energy band, in particular, in the interesting half filling case since it with narrower energy bands. 
Obviously, in a given subspace n, the structure features of energy spectrum of pairing model can be represented by 
the simple representative energy spectrum which consists of a LEL and n(< [L/2]) or L — n(n > [L/2]) RELs. In 
particular, based on the data analysis and numerical fitting for the RELs, we find that the widths of representative 
energy gaps, that is, A\ L \k,n) — E^ L \k, n) — E^ L \k — 1, n) (where eI L \o, 1) = £^(1,1)) from low to high decrease 

near a constant 2V, which means that A^ (ft, n) « Ar (k — l,n) — 2V, and Ar (ft, n) » A^ (ft, n+ 1). From them 
we have 
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Of course , the form of this two additional terms is given by the numerical fitting. Therefore, we obtain the formula 
to estimate the RELs in the subspace n (n < [L/2]) by 
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where A^ L) (1,1) = ^ L) (l,l) -£^(1,1), which can be calculated by the directly diagonlizing the subspace 1 of 
H p . Note that we, based on the data analysis and numerical fitting, add two little addition terms in order to include 
the contribution from high order approximations among the representative energy gaps . When L — 10, A = 20, the 
absolute error of this estimation is ~ 10 -9 and its relative error is ~ 10~ 3 . Then, based on the symmetry in the energy 
spectrum, we can obtain the values of representative energy spectrum in the other subspaces. Fig0| is an example 
which is estimated by our above scheme. The representative energy spectrum in subspaces 11 to 19 has not been 
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FIG. 4: The estimated the lowest energy level and representative energy spectrum (A = 20) for the subspaces from n = 1 to 

71 = 10 

pictured since their shapes can be known by the symmetry of energy spectrum. For the subspace 2, we can calculate 
out the absolute errors of one estimated the lowest energy level and two estimated representative energy levels are 
respectively 5.80231 x 1(T 9 , 2.63859 x 1(T 9 and 6.47206 x 10~ 9 , their relative errors are respectively 0.0000721591, 
0.0000653016 and 0.00147345. Therefore, the lowest energy levels and the representative energy levels can have the 
better and usable precision in numerical calculation as well as estimation by our method. 

It is clear that we have solved the lowest energy levels of all of the subspaces and representative energy spectrum 
of the pairing model in numerical by our direct diagonalization of Fock space. Moreover, we show the structure of 
the energy levels and energy gaps in the pairing model for the reduced model. After finishing the proof of symmetry 
between the energy spectrum of pairing model, we would like to show the origin of leading term of the LEL of H p in 
the near. In principle, our method should be able to extend to some similar spin Hamiltonian systems. Of course, 
because our estimated formula are obtained by the numerical fitting and analytical continuation, it is still a problem 
how large N our scheme is suitable to. In addition, for the width of energy bands we have not found a good estimated 
method yet. More knowledge about the other models needs to explore. 
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